cap program drop adminbunchest
program define adminbunchest, rclass
	args lb threshold
	
	*********************************************************************************
	*	Step 4: Collapse Data into counts of firms by revenue bins					*
	*********************************************************************************
			
					
		collapse (count) adminratio [fw=wt], by(adminscoregroups)   //score is raw score; scoregroups = rounded scores
				
			*** Score polynomial terms 
			
			replace adminscoregroups = 0.1 if adminscoregroups == 0 

			gen adminscore = adminscoregroups
			gen adminscore2 = adminscore^2
			gen adminscore3 = adminscore^3
			gen adminscore4 = adminscore^4
	
			
	*** Estimate Bunching: 
	
	
	mat results = (`lb',0,0,0,0,0,0,0)
	local ub1 = `threshold'+0.2

		forval i = `=`threshold'+.2' (0.2) 18 {
		
			*di "Testing upper bound: `i'"
			
			cap drop exclregn exclregndummy 
			
			qui gen exclregn = adminscore if inrange(adminscoregroups,`lb',`i') // define dummies for bunching region
				qui replace exclregn = 0 if exclregn ==. 
			
			qui egen exclregndummy = group(exclregn)
			qui tab exclregndummy, matrow(regndummynames)
			local maxcats = rowsof(regndummynames) // save dummy name of the last output group
			qui levelsof exclregndummy if exclregn == `threshold', local(notch) //save dummy name of the notch group
			
			*Estimate counterfactual: 
			
			
			qui reg adminratio i.exclregndummy adminscore adminscore2 adminscore3 adminscore4 
			
			* Identify starting and ending columns of results matrix with the coefficient dummies for excluded region
			
			mat A = e(b) // save results matrix
			
			local startcol =  colnumb(A,"2.exclregndummy")
			local notchcol = colnumb(A,"`notch'.exclregndummy")
			local afternotch = `notchcol' + 1
			local endcol =  colnumb(A,"`maxcats'.exclregndummy")   
			
*			di "Start col: `startcol', notchcol: `notchcol', endcol: `endcol'"

			* Add up the missing mass:    
			
			mkmat adminratio if inrange(adminscore,`lb',`=`threshold'+.1')
				mat M_actual = adminratio   // M_actual is 36x1
			
			mkmat adminscore if inrange(adminscore,`lb',`=`threshold'+.1')
				mat Om = adminscore   //Om is 36x1
				
				
			local Mlen = rowsof(Om)   //missing mass length of matrix
			mat C = J(1,`Mlen',1)  // C is 1x36
			
			mata : st_matrix("Om2", st_matrix("Om"):^2)
			mata : st_matrix("Om3", st_matrix("Om"):^3)
			mata : st_matrix("Om4", st_matrix("Om"):^4)

			mat M_est = _b[_cons]*C' + _b[adminscore]*Om + _b[adminscore2]*Om2 + _b[adminscore3]*Om3 + _b[adminscore4]*Om4 
			
				
			mat M_missing = M_est - M_actual // mass in excess of counterfactual density  36x1
			
			mat D = M_missing'*C'  // 1x1
			
			local missingmass = D[1,1]
			 
			* Add up the excess mass:     
			
			mkmat adminratio if inrange(adminscore,`=`threshold'+.1',`i')
				mat B_actual = adminratio    // B_actual is 1x1
			
			mkmat adminscore if inrange(adminscore,`=`threshold'+.1',`i')
				mat Ob = adminscore   // Ob is 1x1
				
						
			local Blen = rowsof(Ob)   //missing mass length of matrix
			mat C = J(1,`Blen',1)    
				
			mata : st_matrix("Ob2", st_matrix("Ob"):^2)
			mata : st_matrix("Ob3", st_matrix("Ob"):^3)
			mata : st_matrix("Ob4", st_matrix("Ob"):^4)

			mat B_est = _b[_cons]*C' + _b[adminscore]*Ob + _b[adminscore2]*Ob2 + _b[adminscore3]*Ob3 + _b[adminscore4]*Ob4 
			
			
			mat B_excess = B_actual - B_est  // mass in excess of counterfactual density
			
			mat D = B_excess'*C'
			
			local excessmass = D[1,1]
			
			* Difference between missing and excess mass: 
			
			local diff = `excessmass' - `missingmass'
			
			di "Difference between excess and missing mass is: `diff' for upper bound `i'"
			
			* Counterfactual densities at the notch and at the upper bound:
				
			local cflb = _b[_cons] + `lb'*_b[adminscore]+ (`lb'^2)*_b[adminscore2] + (`lb'^3)* _b[adminscore3] + (`lb'^4)* _b[adminscore4] 
				
			local cfnotch = _b[_cons] + `threshold'*_b[adminscore]+ (`threshold'^2)*_b[adminscore2] + (`threshold'^3)* _b[adminscore3] + (`threshold'^4)* _b[adminscore4] 
				
			local cfub = _b[_cons] + `i'*_b[adminscore]+ (`i'^2)*_b[adminscore2] + (`i'^3)* _b[adminscore3] + (`i'^4)* _b[adminscore4] 

			* Estimated average bunching: 
			
			local b_av = `excessmass'/(0.5*(`cflb'+`cfub'))
			
			mat resultsappend = (`lb',`i',`excessmass',`missingmass',`diff',`cflb',`cfub',`b_av')
			mat results = results \ resultsappend
			
		}
		
		mat colnames results = lb ub excess missing diff cflb cfub b_av
		
		clear 
		svmat results, names(col)
		
		* Keep the bunching estimate for the lowest difference between excess and missing mass: 
		
		drop if ub == 0 
		gen absdiff = abs(diff)
		sort absdiff
		
		keep if _n == 1
		mkmat _all, matrix(R)
				
		
		return scalar ub_est = R[1,2]
		return scalar excess_est = R[1,3]
		return scalar b_av_est = R[1,8]  
end
